%% Carbon Taxes Around the World: Cooperation, Strategic Interactions, and Spillovers
% IMF Economic Review
% Alessandro Moro and Valerio Nispi Landi
% Replication files
% This file simulates baseline model

close all;

warning off
LIN=0;         % if 1 the tax increases linearly over time
DAM=0;         % if 1 damage is higher
LOOP=0;        % 1=run the digamma loop; 0=run a single simulation

%%%%%%%%%%%%%%%%%%%%%%
T=401;              % periods of simulation
TB=180;             % from this period on, all emissions are abated in the baseline police
TO=380;             % from this period on, all emissions are abated in the other policy

if LOOP==1
DG= 0.0000:0.0001:0.007;   
DGZ=0.0000:0.0001:0.007; 
nn=length(DG);
nnz=length(DGZ);
DGend=0.02;
DGZend=0.02;
DG(nn+1:nn+1000*(DGend-DG(end)))=DG(end)+0.001:0.001:DGend;
DGZ(nnz+1:nnz+1000*(DGZend-DGZ(end)))=DGZ(end)+0.001:0.001:DGZend;
%DGZ=0; 

else
load digs    
DG= digamma;   
DGZ=digammaz;  
end

if LOOP==1 & LIN==1
DG=DG/10;
DGZ=DGZ/10;
elseif LOOP==1 & DAM==1
DG=10*DG;
DGZ=10*DGZ;
end

% This matrix gives a welfare level for every combination of DG-DGZ. It will be filled by the loop
WW= -1e8*ones(length(DG),length(DGZ));
WWz=-1e8*ones(length(DG),length(DGZ));
WWw=-1e8*ones(length(DG),length(DGZ));

N=(length(DG)*length(DGZ));

cc=0;                % Counter

for kk=1:length(DG)   % tau loop starts
    digamma=DG(kk);
for mm=1:length(DGZ)  % tauz loop starts
    digammaz=DGZ(mm);


tt=zeros(T,1);
ttz=zeros(T,1);
tt(1)=tauss;
ttz(1)=tauzss;
tt(TO+1:end)=tau_clean;     % from 2201 on all emissions are abated
ttz(TO+1:end)=tauz_clean;   % from 2201 on all emissions are abated
zs=(iota.^(0:1:T-2))';

%BASE;         % 1=baseline policy; 0=other (optimal) policy, 
%2: only country F follows the baseline policy; 3: both countries set the tax always equal to the ss

% COUNTRY 1 BASELINE, COUNTRY 2 BASELINE
if BASE==1
 for j=2:TB
 tt(j)=1.01*tt(j-1); 
 ttz(j)=1.01*ttz(j-1);% before TB price of carbon grows at 1%
 end
tt(TB+1:end)=tau_clean; 
ttz(TB+1:end)=tauz_clean; 

% COUNTRY 1 OPTIMAL, COUNTRY 2 BASELINE
 elseif BASE==2  % only country 2 follows the baseline policy
 for j=2:TB
 ttz(j)=1.01*ttz(j-1);% before TB price of carbon grows at 1%
 end
ttz(TB+1:end)=tauz_clean; 

   if LIN==0 
   tt(1:TO)=[tauss; min(tau_clean,tauss+digamma*zs(1:TO-1))];
   elseif LIN==1
   tt(1:TO)=[tauss; min(tau_clean,tauss+digamma*(1:1:TO-1)')];
   end

% COUNTRY 1 CONSTANT, COUNTRY 2 CONSTANT
elseif BASE==3
tt(1:TO)=tauss*ones(TO,1);
ttz(1:TO)=tauss*ones(TO,1);

% BASE=0: COUNTRY 1 OPTIMAL, COUNTRY 2 OPTIMAL

else
if LIN==0
% First value of the tax is the steady state, then the tax is tauss+digamma*z0, tauss+ digamma*z1, etc
tt(1:TO)=[tauss; min(tau_clean,tauss+digamma*zs(1:TO-1))];
ttz(1:TO)=[tauzss;min(tauz_clean,tauzss+digammaz*zs(1:TO-1))];
elseif LIN==1;
tt(1:TO)=[tauss; min(tau_clean,tauss+digamma*(1:1:TO-1)')];
ttz(1:TO)=[tauzss;min(tauz_clean,tauzss+digammaz*(1:1:TO-1)')];  
end
end

%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%Endogenous Variables %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
var
%AE     %EME
c       cz          % consumption
yH      yFz         % production
i       i_zz        % investment (domestic)
i_zu    i_zd        % investment (foreign)
kG      kG_zz       % green domestic capital 
kB      kB_zz       % brown foreign capital 
kG_zu   kG_zd       % green foreign capital
kB_zu   kB_zd       % brown foreign capital 
yB      yBz         % brown production
yG      yGz         % green production
hB      hBz         % brown hours
hG      hGz         % green hours
h       hz          % total hours
e       ez          % emissions
mu      muz         % abatement
w       wz          % wage
r       rz          % risk-free rate
rkG     rkG_zz      % rental rate green domestic capital
rkB     rkB_zz      % rental rate brown domestic capital
rkG_zu  rkG_zd      % rental rate brown foreign capital
rkB_zu  rkB_zd      % rental rate green domestic capital
ErkG     ErkG_zz    % Expected rental rate green domestic capital
ErkB     ErkB_zz    % Expected rental rate brown domestic capital
ErkG_zu  ErkG_zd    % Expected rental rate brown foreign capital
ErkB_zu  ErkB_zd    % Expected rental rate green domestic capital

pB      pBz         % brown price
pG      pGz         % green price
b                   % financial position of households
dF                  % financial position of EME banks
dH                  % financial position of AE banks
pH                  % price of AE good in terms of AE CPI
pF                  % price of EME good in terms of AE CPI
s                   % real exchange rate (if it goes up, AE currency depreciates)
Ab     Abz          % abatement spending
pB2    pB2z         % net brown price
gdp    gdpz         % GDP
mp     mpz          % imports
xp     xpz          % exports    
tb     tbz          % trade balance
uip          
proof
proof2
proof3
;

%%
%%%%%%%%%%%%%%%%%%%%%%%Exogenous Variables%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
varexo 
tau    tauz      % tax shock

;   
    
%%
%%%%%%%%%%%%%%%%%%%%%%%Parameters%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

parameters
gamma
gammaz
eta
chi
beta
betaz
GAM
GAMz
n
delta
deltax
iota
zeta
csi
kappaA
kappaAz
nu
v
g
gz
alpha1
alpha2
alpha1z
alpha2z
A
Az
yHss
yFzss
pBss
pBzss
pHss
bss
kBss
k_zdss
k_zuss
rkGss
rkG_zzss
tauss
tauzss
iss
i_zzss
yHend 
pHend 
pBend 
pBzend
yFzend
bend 
kBend
rkGend 
rkG_zzend
iend 
i_zzend
cend 
czend 
kappaX
kappaXz
psiX
tauend
tauzend
zetae 
zetaez
x1_base
;

load par;    % load mat file created in console
load paropt; % load mat file created in console_final

for jj=1:length(M_.param_names)
set_param_value(M_.param_names{jj},eval(M_.param_names{jj})); 
end;



%%
%%%%%%%%%%%%%%%%%%%%%%%Non-Linear Model%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

model; 
% Households
1=c/c(1) *beta  /iota*r;                                                   %(1)
1=cz/cz(1)*betaz/iota*rz;                                                  %(2)
1=beta /iota*c/c(1)*(rkG(1)   +1-delta);                                   %(3)
1=betaz/iota*cz/cz(1)*(rkG_zz(1)+1-delta);                                 %(4)
1=beta /iota*c/c(1) *(rkB(1)  +1-delta);                                   %(5)
1=betaz/iota*cz/cz(1)*(rkB_zz(1)+1-delta);                                 %(6)
1=beta /iota*c/c(1)*s(1)/s*(rkG_zd(1)+1-delta);                            %(7)
1=betaz/iota*cz/cz(1)*s/s(1)*(rkG_zu(1)+1-delta);                          %(8)
1=beta /iota*c/c(1)  *s(1)/s*(rkB_zd(1)+1-delta);                          %(9)
1=betaz/iota*cz/cz(1)*s/s(1)*(rkB_zu(1)+1-delta);                          %(10)
kG+kB=(1-delta)*(kG(-1)+kB(-1))/iota+i;                                    %(11)                  
kG_zz+kB_zz=(1-delta)*(kG_zz(-1)+kB_zz(-1))/iota+i_zz;                     %(12)                                                   
kG_zd+kB_zd=(1-delta)*(kG_zd(-1)+kB_zd(-1))/iota+i_zd;                     %(13)                                  
kG_zu+kB_zu=(1-delta)*(kG_zu(-1)+kB_zu(-1))/iota+i_zu;                     %(14)   
c*h^(chi)=w;                                                               %(15)
cz*hz^(chi)=wz;                                                            %(16)
% Brown firms
yB= A* ((kB(-1)   /iota)^alpha1) *((kB_zu(-1)/iota)^alpha2) *(hB^ (1-alpha1-alpha2));   %(17)
yBz=Az*((kB_zz(-1)/iota)^alpha1z)*((kB_zd(-1)/iota)^alpha2z)*(hBz^(1-alpha1z-alpha2z)); %(18)
alpha1* pB2 *yB =kB(-1)   /iota*rkB;                                      %(19)
alpha1z*pB2z*yBz=kB_zz(-1)/iota*rkB_zz;                                    %(20)
alpha2* pB2* yB =kB_zu(-1)/iota*rkB_zu;                                    %(21)
alpha2z*pB2z*yBz=kB_zd(-1)/iota*rkB_zd;                                    %(22)
(1-alpha1-alpha2)  *pB2 *yB =hB*w;                                         %(23)
(1-alpha1z-alpha2z)*pB2z*yBz=hBz*wz;                                       %(24)
mu^nu= min(1,(tau*zetae  /kappaA));                                        %(25)
muz^nu=min(1,(tauz*zetaez/kappaAz));                                       %(26)
% Green firms
yG =A* ((kG(-1)   /iota)^alpha1) *((kG_zu(-1)/iota)^alpha2) *(hG ^(1-alpha1-alpha2));   %(27)
yGz=Az*((kG_zz(-1)/iota)^alpha1z)*((kG_zd(-1)/iota)^alpha2z)*(hGz^(1-alpha1z-alpha2z)); %(28)
alpha1* pG* yG =kG(-1)   /iota*rkG;                                        %(29)
alpha1z*pGz*yGz=kG_zz(-1)/iota*rkG_zz;                                     %(30)
alpha2* pG* yG =kG_zu(-1)/iota*rkG_zu;                                     %(31)
alpha2z*pGz*yGz=kG_zd(-1)/iota*rkG_zd;                                     %(32)
(1-alpha1-alpha2)  *pG*yG  =hG*w;                                          %(33)
(1-alpha1z-alpha2z)*pGz*yGz=hGz*wz;                                        %(34)
% Final-good firms
yB =zeta*(pB/pH)   ^(-csi)*yH;                                             %(35)
yBz=zeta*(s*pBz/pF)^(-csi)*yFz;                                            %(36)
yG= (1-zeta)*(pG/pH)   ^(-csi)*yH;                                         %(37)
yGz=(1-zeta)*(s*pGz/pF)^(-csi)*yFz;                                        %(38)
pH*yH   =pG*yG+pB*yB;                                                      %(39)
pF*yFz/s=pGz*yGz+pBz*yBz;                                                  %(40)
1=(1-gamma) *pH^(1-eta)+gamma *pF^(1-eta);                                 %(41)
s^(1-eta)=(1-gammaz)*pF^(1-eta)+gammaz*pH^(1-eta);                         %(42)
% Banks
b=-(1/GAM*         (beta/iota*c/c(1)  *(r-rz*s(1)/s))+(1-n)/n*dF+v);       %(43)
b= (s/GAMz*(1-n)/n*(betaz/iota*cz/cz(1)*(rz-r*s/s(1)))-(dH+v));            %(44)
% Market clearing
0=(1-n)/n*dF+dH+v+b;                                                       %(45)
h= hG +hB;                                                                 %(46)
hz=hGz+hBz;                                                                %(47)
n*yH     =n*     ((1-gamma)* pH   ^(-eta)*(c+ i   +i_zu+Ab) +g) 
+(1-n)*(gammaz*(pH/s)^(-eta)*(cz+i_zz+i_zd+Abz));                          %(48)
(1-n)*yFz=(1-n)*((1-gammaz)*(pF/s)^(-eta)*(cz+i_zz+i_zd+Abz)+gz)
+n    *gamma  *pF^(-eta)    *(c  +i  +i_zu+Ab);                            %(49)
pH*yH=c+i+pH*g+Ab+b+(1-n)/n*s*i_zd-s*(1-n)/n*(rkB_zd*kB_zd(-1)+rkG_zd*kG_zd(-1))/iota
+rz(-1)*s/s(-1)*(dH(-1)+v)/iota+(rkB_zu*kB_zu(-1)
+rkG_zu*kG_zu(-1))/iota-r(-1)/iota*(b(-1)+dH(-1)+v);                       %(50)
% Emissions
e= zetae*(1-mu)*yB;                                                        %(51)
ez=zetaez*(1-muz)*yBz;                                                     %(52)

% Auxiliary variables
Ab=kappaA /(1+nu)*(mu ^(1+nu))*yB;                                         
Abz=kappaAz/(1+nu)*(muz^(1+nu))*yBz;                                       
pB2 =pB -tau *(1-mu)*zetae -kappaA/(1+nu)*mu ^(1+nu);                      
pB2z=pBz-tauz*(1-muz)*zetaez-kappaAz/(1+nu)*muz^(1+nu);   
gdp=pH*yH;
gdpz=pF*yFz/s;
mp =pF*gamma   *pF    ^(-eta)*(c +i   +i_zu+Ab);
mpz=pH/s*gammaz*(pH/s)^(-eta)*(cz+i_zz+i_zd+Abz);
xp=pH*gammaz*(pH/s)^(-eta)*(1-n)/n*(cz+i_zz+i_zd+Abz);
xpz=n/(1-n)*pF/s*gamma *pF ^(-eta)*(c+i +i_zu+Ab);
tb =xp -mp;
tbz=xpz-mpz;
uip=rz*s(1)/s-r;
ErkG=rkG(1);
ErkG_zz=rkG_zz(1);
ErkB=rkB(1);
ErkB_zz=rkB_zz(1);
ErkG_zu=rkG_zu(1);
ErkG_zd=rkG_zd(1);
ErkB_zu=rkB_zu(1);
ErkB_zd=rkB_zd(1);  
proof=cz+i_zd+i_zz+pF/s*gz+kappaAz/(1+nu)*muz^(1+nu)*yBz+tbz-(pF/s*yFz);
proof2=pF/s*yFz-(cz+i_zz+pF*gz/s+Abz
       -n/((1-n)*s)*b+n/(s*(1-n))*i_zu
       +rkB_zd*kB_zd(-1)/iota+rkG_zd*kG_zd(-1)/iota+
       -n/((1-n)*s)*(rkB_zu*kB_zu(-1)/iota+rkG_zu*kG_zu(-1)/iota)
       +r(-1)/s*n/(1-n)*(b(-1)+dH(-1)+v)/iota-n/(1-n)*rz(-1)/s(-1)*(dH(-1)+v)/iota);
proof3=gdp-(c+i+pH*g+Ab
       +b+(1-n)/n*s*i_zd
       -s*(1-n)/n*(rkB_zd*kB_zd(-1)/iota+rkG_zd*kG_zd(-1)/iota)
       +(rkB_zu*kB_zu(-1)/iota+rkG_zu*kG_zu(-1)/iota)
       -r(-1)*(b(-1)+dH(-1)+v)/iota+rz(-1)*s/s(-1)*(dH(-1)+v)/iota);       
end;


initval;
yH=yHss;
yFz=yFzss;
pB=pBss;
pBz=pBzss;
pH=pHss;
b=bss;
kB=kBss;
tau=tauss;
tauz=tauzss;
mu= min(1,(tau*zetae /kappaA)^(1/nu));                                                  
muz=min(1,(tauz*zetaez/kappaAz)^(1/nu));                                                       
i=iss;
i_zz=i_zzss;
rkB=iota/beta-(1-delta);            
rkB_zd=iota/beta-(1-delta);         
rkB_zu=iota/betaz-(1-delta);       
rkB_zz=iota/betaz-(1-delta);       
rkG=iota/beta-(1-delta);           
rkG_zz=iota/betaz-(1-delta);       
rkG_zd=iota/beta-(1-delta);         
rkG_zu=iota/betaz-(1-delta);      
r=iota/beta;
rz=iota/betaz;
pF=((1-(1-gamma)*pH^(1-eta))/gamma)^(1/(1-eta));
s=(gammaz*pH^(1-eta)+(1-gammaz)*pF^(1-eta))^(1/(1-eta));
dH=(1-beta/betaz)/GAM;
dF=s/GAMz*(betaz/beta-1);
i_zd=k_zdss*(1-(1-delta)/iota);
i_zu=k_zuss*(1-(1-delta)/iota);
yB=zeta*(pB/pH)^(-csi)*yH;
pB2=pB-tau*(1-mu)*zetae-kappaA/(1+nu)*mu^(1+nu);
pB2z=pBz-tauz*(1-muz)*zetaez-kappaAz/(1+nu)*muz^(1+nu);
kG=alpha1*iota/rkG*(pH*yH-pB*yB);
kB_zu=alpha2*iota/rkB_zu*pB2*yB;
kG_zu=alpha2*iota/rkG_zu*(pH*yH-pB*yB);
hB=(yB/(A*(kB/iota)^alpha1*(kB_zu/iota)^alpha2))^(1/(1-alpha1-alpha2));
w=(1-alpha1-alpha2)*pB2*yB/hB;
hG=(1-alpha1-alpha2)*(pH*yH-pB*yB)/w;
h=hG+hB;
yG=A*(kG/iota)^alpha1*(kG_zu/iota)^alpha2*(hG)^(1-alpha1-alpha2);
pG=(pH*yH-pB*yB)/yG;
yBz=zeta*(s*pBz/pF)^(-csi)*yFz;
kB_zz=alpha1z*iota*pB2z*yBz/rkB_zz;
kG_zz=alpha1z*iota*(pF*yFz/s-pBz*yBz)/rkG_zz;
kB_zd=alpha2z*iota*pB2z*yBz/rkB_zd;
kG_zd=alpha2z*iota*(pF*yFz/s-pBz*yBz)/rkG_zd;
gdp=pH*yH;
c=-(i+pH*g+kappaA/(1+nu)*mu^(1+nu)*yB+b+(1-n)/n*s*i_zd-s*(1-n)/n*(rkB_zd*kB_zd    
+rkG_zd*kG_zd)    /iota+(rkB_zu*kB_zu    +rkG_zu*kG_zu)    /iota
+rz            *(dH+v)    /iota-r/iota    *(b+    dH    +v)-gdp);
cz=n/((1-n)*gammaz)*(pH/s)^eta*(yH-((1-gamma)*pH^(-eta)*(c+i+kappaA/(1+nu)*mu^(1+nu)*yB+i_zu)+g))
-(i_zz+i_zd+kappaAz/(1+nu)*muz^(1+nu)*yBz);
wz=(cz^(1/chi)*(1-alpha1z-alpha2z)*(pB2z*yBz+pGYz))^(chi/(1+chi));
hz=(wz/cz)^(1/chi);
hBz=(1-alpha1z-alpha2z)*pB2z*yBz/wz;
hGz=hz-hBz;
yGz=Az*(kG_zz/iota)^alpha1z*(kG_zd/iota)^alpha2z*(hGz)^(1-alpha1z-alpha2z);
pGz=(pF*yFz/s-pBz*yBz)/yGz;
e=zetae*(1-mu)*yB;                                                                                                 
ez=zetaez*(1-muz)*yBz;   
Ab=kappaA/(1+nu)*(mu^(1+nu))*yB;                                                                                           
Abz=kappaAz/(1+nu)*(muz^(1+nu))*yBz; 
gdpz=pF*yFz/s;
mp =pF*gamma   *pF    ^(-eta)*(c +i   +i_zu+Ab);
mpz=pH/s*gammaz*(pH/s)^(-eta)*(cz+i_zz+i_zd+Abz);
xp=pH*gammaz*(pH/s)^(-eta)*(1-n)/n*(cz+i_zz+i_zd+Abz);
xpz=n/(1-n)*pF/s*gamma *pF ^(-eta)*(c+i +i_zu+Ab);
tb =xp -mp;
tbz=xpz-mpz; 
uip=rz-r;
ErkG=rkG;
ErkG_zz=rkG_zz;
ErkB=rkB;
ErkB_zz=rkB_zz;
ErkG_zu=rkG_zu;
ErkG_zd=rkG_zd;
ErkB_zu=rkB_zu;
ErkB_zd=rkB_zd; 
proof2=0;
proof3=0;
end;

%steady;

endval;
yH=yHend;
yFz=yFzend;
pB=pBend;
pBz=pBzend;
pH=pHend;
b=bend;
kB=kBend;
r=iota/beta;
rz=iota/betaz;         
rkB=iota/beta-(1-delta);            
rkB_zd=iota/beta-(1-delta);        
rkB_zu=iota/betaz-(1-delta);     
rkB_zz=iota/betaz-(1-delta);       
rkG=iota/beta-(1-delta);           
rkG_zz=iota/betaz-(1-delta);       
rkG_zd=iota/beta-(1-delta);         
rkG_zu=iota/betaz-(1-delta);       
tau=tauend;
tauz=tauzend;
mu= min(1,(tau*zetae /kappaA)^(1/nu));                                                  
muz=min(1,(tauz*zetaez/kappaAz)^(1/nu));                                                     
i=iend;
i_zz=i_zzend;
pF=((1-(1-gamma)*pH^(1-eta))/gamma)^(1/(1-eta));
s=(gammaz*pH^(1-eta)+(1-gammaz)*pF^(1-eta))^(1/(1-eta));
dH=(1-beta/betaz)/GAM;
dF=s/GAMz*(betaz/beta-1);
yB=zeta*(pB/pH)^(-csi)*yH;
pB2=pB-tau*(1-mu)*zetae-kappaA/(1+nu)*mu^(1+nu);
pB2z=pBz-tauz*(1-muz)*zetaez-kappaAz/(1+nu)*muz^(1+nu);
kG=alpha1*iota/rkG*(pH*yH-pB*yB);
kB_zu=alpha2*iota/rkB_zu*pB2*yB;
kG_zu=alpha2*iota/rkG_zu*(pH*yH-pB*yB);
i_zu=(kB_zu+kG_zu)*(1-(1-delta)/iota);
hB=(yB/(A*(kB/iota)^alpha1*(kB_zu/iota)^alpha2))^(1/(1-alpha1-alpha2));
w=(1-alpha1-alpha2)*pB2*yB/hB;
hG=(1-alpha1-alpha2)*(pH*yH-pB*yB)/w;
h=hG+hB;
yG=A*(kG/iota)^alpha1*(kG_zu/iota)^alpha2*(hG)^(1-alpha1-alpha2);
pG=(pH*yH-pB*yB)/yG;
yBz=zeta*(s*pBz/pF)^(-csi)*yFz;
gdpz=pF*yFz/s;
kB_zz=alpha1z*iota*pB2z*yBz/rkB_zz;
kG_zz=alpha1z*iota*(gdpz-pBz*yBz)/rkG_zz;
kB_zd=alpha2z*iota*pB2z*yBz/rkB_zd;
kG_zd=alpha2z*iota*(gdpz-pBz*yBz)/rkG_zd;
i_zd=(kB_zd+kG_zd)*(1-(1-delta)/iota);
gdp=pH*yH;
c=-(i+pH*g+kappaA/(1+nu)*mu^(1+nu)*yB+b+(1-n)/n*s*i_zd-s*(1-n)/n*(rkB_zd*kB_zd    
+rkG_zd*kG_zd)    /iota+(rkB_zu*kB_zu    +rkG_zu*kG_zu)    /iota
+rz            *(dH+v)    /iota-r/iota    *(b+    dH    +v)-gdp);
cz=n/((1-n)*gammaz)*(pH/s)^(eta)*(yH-((1-gamma)*pH^(-eta)*(c+i+kappaA/(1+nu)*mu^(1+nu)*yB+i_zu)+g))
-(i_zz+i_zd+kappaAz/(1+nu)*muz^(1+nu)*yBz);
wz=(cz^(1/chi)*(1-alpha1z-alpha2z)*(pB2z*yBz+pGYz))^(chi/(1+chi));
hz=(wz/cz)^(1/chi);
hBz=(1-alpha1z-alpha2z)*pB2z*yBz/wz;
hGz=hz-hBz;
yGz=Az*(kG_zz/iota)^alpha1z*(kG_zd/iota)^alpha2z*(hGz)^(1-alpha1z-alpha2z);
pGz=(pF*yFz/s-pBz*yBz)/yGz;
e=zetae*(1-mu)*yB;                                                                                                 
ez=zetaez*(1-muz)*yBz;   
Ab=kappaA/(1+nu)*(mu^(1+nu))*yB;                                                                                           
Abz=kappaAz/(1+nu)*(muz^(1+nu))*yBz; 
mp =pF*gamma   *pF    ^(-eta)*(c +i   +i_zu+Ab);
mpz=pH/s*gammaz*(pH/s)^(-eta)*(cz+i_zz+i_zd+Abz);
xp=pH*gammaz*(pH/s)^(-eta)*(1-n)/n*(cz+i_zz+i_zd+Abz);
xpz=n/(1-n)*pF/s*gamma *pF ^(-eta)*(c+i +i_zu+Ab);
tb =xp -mp;
tbz=xpz-mpz;
uip=rz-r;
ErkG=rkG;
ErkG_zz=rkG_zz;
ErkB=rkB;
ErkB_zz=rkB_zz;
ErkG_zu=rkG_zu;
ErkG_zd=rkG_zd;
ErkB_zu=rkB_zu;
ErkB_zd=rkB_zd; 
proof2=0;
proof3=0;
end;
%steady;


% Start from period 2, the first period is the SS and is added by Dynare
t2=tt(2:end);
t2z=ttz(2:end);


shocks;
var tau;
periods 1:400;
values (t2);
var tauz;
periods 1:400;
values (t2z);
end;


%% Simulation
perfect_foresight_setup(periods=400);
perfect_foresight_solver(stack_solve_algo=7, solve_algo=10); 



if oo_.deterministic_simulation.status ~= 1
disp(['Computation fails in Dynare']);
return
end

%% Verify steady state
ver1=abs(c(1)-css);
ver2=abs(yFz(end)-yFzend);

if ver1>1e-09
disp(['*****!!!!!!Initial steady states do not coincide'])
return
ERROR
elseif ver2>1e-06
disp(['*****!!!!!!Final steady states do not coincide'])
return
ERROR
end    

%% Include productivity
zs=(iota.^(0:1:T))';
cs=c.*zs;
czs=cz.*zs;
es=e.*zs;
ezs=ez.*zs;
yBzs=yBz.*zs;
gdps=gdp.*zs;
gdpzs=gdpz.*zs;
kBs=kB.*zs;
kGs=kG.*zs;
kB_zus=kB_zu.*zs;
kG_zus=kG_zu.*zs;
kB_zzs=kB_zz.*zs;
kG_zzs=kG_zz.*zs;
kB_zds=kB_zd.*zs;
kG_zds=kG_zd.*zs;
ErkBs=ErkB;
ErkGs=ErkG;
ErkB_zus=ErkB_zu;
ErkG_zus=ErkG_zu;
ErkB_zzs=ErkB_zz;
ErkG_zzs=ErkG_zz;
ErkB_zds=ErkB_zd;
ErkG_zds=ErkG_zd;
ss=s;
bYs=b./gdp;
tbYs=tb./gdp;
eWs=n*es+(1-n)*ezs;
rs=r;
rzs=rz;
mus=mu;
muzs=muz;
uips=uip;
prices= oo_.exo_simul(:,1)*conv_gdp/3.67;
pricezs=oo_.exo_simul(:,2).*ss*conv_gdp/3.67;
priceWs=n*prices+(1-n)*pricezs;
muW=n*mu+(1-n)*muz;
Etot=n*es+(1-n)*ezs;
x1s=zeros(402,1);
x2s=zeros(402,1);
x1s(1)=x1_base;             
x2s(1)=xGtC-x1_base;
for j=2:402
 x1s(j)=x1s(j-1)+phiL*Etot(j);
 x2s(j)=(1-phi)*x2s(j-1)+phi0*(1-phiL)*Etot(j);
end
xs=x1s+x2s;

%% Welfare
temps=3*log(xs./xbar)./log(2);
BETA=(beta.^(0:1:T))';
BETAz=(betaz.^(0:1:T))';

Us= BETA.*(log(cs)-1/(1+chi)*h.^(1+chi)-(kappaX/(1+psiX)*(xs-xbar).^(1+psiX)));
Uzs=BETAz.*(log(czs)-1/(1+chi)*hz.^(1+chi)-(kappaXz/(1+psiX)*(xs-xbar).^(1+psiX)));


W=sum(Us(2:T)); 
Wz=sum(Uzs(2:T)); 
Ww=n*W+(1-n)*Wz;

%%%

WW(kk,mm)=W;
WWz(kk,mm)=Wz;
WWw(kk,mm)=Ww;   
  
cc=cc+1;        % update the counter
disp(['-----------------------------------------------------------------------']);
disp(['Percentage completed: ',num2str(100*cc/N),'%']);
disp(['-----------------------------------------------------------------------']);



end   % tau loop ends
end   % tauz  loop ends


if max(DGZ)==0
    [~, ind_max] = max(WW);
    DG_opt=DG(ind_max);
    fprintf('Optimal digamma when EME follows baseline policy: %.4f\n', DG_opt);
end